function [act logact] = kmeansactivRobust(mix, data)


K = mix.ncentres;
[N D] = size(data.Y);

% precomputed values and vectorised data
Dlog2pi = D*log(2*pi);
invC_vec = reshape(data.invC, D*D, N);
logdetC = data.logdetC;
invCy = data.invCy;
yinvCy = data.yinvCy;

% vectorised model params
mu = mix.centres;
mumut = zeros(D, D, K);
for k = 1:K
    mumut(:, :, k) = mu(k, :)'*mu(k, :);
end
mumut_vec = reshape(mumut, D*D, K);

% compute log act using vectorised forms
logact = -Dlog2pi - logdetC - yinvCy;   % yinvCy is part of the quadratic form
logact = logact*ones(1, K);
logact = logact - invC_vec'*mumut_vec + 2*invCy*mu';
logact = 0.5*logact;

act = exp(logact);
